El ecosistema geoespacial de R

Trabajo previo

Lea los capítulos del 1 al 2 de Lovelace, R., Nowosad, J. & Muenchow, J. (2020). Geocomputation with R.

Preparativos

Carga de paquetes

# Paquete para manipulación de datos
library(dplyr)

# Paquete para manejo de datos vectoriales
library(sf)

# Paquete para manejo de datos raster
library(terra)

# Paquete para mapas interactivos
library(leaflet)

Conjuntos de datos

Introducción

La comunidad de programadores de R ha desarrollado un conjunto de paquetes para el manejo de datos geoespaciales, tanto en formatos vectoriales como raster. Algunos de los principales de estos paquetes son:

Algunos paquetes de graficación, como ggplot2, también cuentan con algunas capacidades para procesamiento de datos geoespaciales.

En CRAN Task View: Analysis of Spatial Data, puede encontrarse un resumen detallado de los paquetes geoespaciales de R.

Datos vectoriales

El modelo vectorial

El modelo vectorial de datos está basado en puntos localizados en un sistema de referencia de coordenadas (CRS). Los puntos individuales pueden representar objetos independientes (ej. la localización de un poste eléctrico o de una cabina telefónica) o pueden también agruparse para formar geometrías más complejas como líneas o polígonos. Por lo general, los puntos tienen solo dos dimensiones (x, y), a las que se les puede agregar una tercera dimensión z, usualmente correspondiente a la altitud sobre el nivel del mar.

El estándar Simple Features

Simple Features (o Simple Feature Access) es un estándar abierto de la Organización Internacional de Estandarización (ISO) y del Open Geospatial Consortium (OGC) que especifica un modelo común de almacenamiento y acceso para geometrías de dos dimensiones (líneas, polígonos, multilíneas, multipolígonos, etc.). El estándar es implementado por muchas bibliotecas y bases de datos geoespaciales como sf, GDAL, PostgreSQL/PostGIS, SQLite/SpatiaLite, Oracle Spatial y Microsoft SQL Server, entre muchas otras.

La especificación define 17 tipos de geometrías, de las cuales siete son las más comúnmente utilizadas. Estas últimas se muestran en la figura 1.

Figura 1. Tipos de geometrías de Simple Features más usadas. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#vector-data)

El paquete sf

El paquete sf (de Simple Features) de R implementa los modelos de datos de las geometrías de tipo vectorial: puntos, líneas, polígonos, sus versiones múltiples y las colecciones de geometrías. Está basado en bibliotecas de sofware ampliamente utilizadas en aplicaciones geoespaciales:

sf provee acceso, desde un mismo paquete de R, a la funcionalidad de estas tres bibliotecas, proporcionando así una interfaz unificada para leer y escribir datos geoespaciales mediante GDAL, realizar operaciones con geometrías mediante GEOS y efectuar transformaciones entre sistemas de coordenadas mediante PROJ.

En sf, los conjuntos de datos geoespaciales se almacenan en objetos de una clase también llamada sf, los cuales son data frames que contiene una columna especial para las geometrías. Esta columna se denomina generalmente geom o geometry (aunque pueden tener cualquier otro nombre). El manejo de datos geoespaciales como data frames permite manipularlos con las funciones ya desarrolladas para este tipo de datos y con la misma forma de referenciar las filas (observaciones) y las columnas (variables).

Métodos del paquete sf

La lista de métodos de la clase sf puede obtenerse a través de la función methods():

# Métodos de la clase sf
methods(class = "sf")
 [1] [                     [[<-                  $<-                  
 [4] aggregate             anti_join             arrange              
 [7] as.data.frame         cbind                 coerce               
[10] dbDataType            dbWriteTable          distinct             
[13] dplyr_reconstruct     ext                   extract              
[16] filter                full_join             group_by             
[19] group_split           identify              initialize           
[22] inner_join            left_join             mask                 
[25] merge                 mutate                plot                 
[28] print                 rasterize             rbind                
[31] rename                right_join            rowwise              
[34] sample_frac           sample_n              select               
[37] semi_join             show                  slice                
[40] slotsFromS3           st_agr                st_agr<-             
[43] st_area               st_as_s2              st_as_sf             
[46] st_bbox               st_boundary           st_buffer            
[49] st_cast               st_centroid           st_collection_extract
[52] st_convex_hull        st_coordinates        st_crop              
[55] st_crs                st_crs<-              st_difference        
[58] st_filter             st_geometry           st_geometry<-        
[61] st_inscribed_circle   st_interpolate_aw     st_intersection      
[64] st_intersects         st_is_valid           st_is                
[67] st_join               st_line_merge         st_m_range           
[70] st_make_valid         st_nearest_points     st_node              
[73] st_normalize          st_point_on_surface   st_polygonize        
[76] st_precision          st_reverse            st_sample            
[79] st_segmentize         st_set_precision      st_shift_longitude   
[82] st_simplify           st_snap               st_sym_difference    
[85] st_transform          st_triangulate        st_union             
[88] st_voronoi            st_wrap_dateline      st_write             
[91] st_z_range            st_zm                 summarise            
[94] transform             transmute             ungroup              
[97] vect                 
see '?methods' for accessing help and source code

Seguidamente, se describen y ejemplifican algunos de los métodos básicos de la clase sf.

st_read() - lectura de datos

La función st_read() lee datos vectoriales de una fuente en formato geoespacial (ej. shapefiles, archivos GeoJSON, bases de datos geoespaciales) y los recupera en un objeto sf.

# Lectura de una capa vectorial (GeoJSON) de provincias de Costa Rica
provincias <-
  st_read(
    "https://github.com/tpb728O-programaciongeoespacialr/2021ii/raw/main/datos/ign/delimitacion-territorial-administrativa/provincias.geojson",
    quiet = TRUE
  )

st_read() también puede crear objetos sf a partir de archivos de texto.

# Lectura de un archivo CSV con registros de presencia de felinos en Costa Rica
felidae <-
  st_read(
    "/home/mfvargas/tpb728O-programaciongeoespacialr/2021ii/datos/gbif/felidae.csv",
    drivers = "CSV", # esta opción parece ser necesario en las últimas versiones de GDAL cuando se lee a través de HTTP
    options = c(
      "X_POSSIBLE_NAMES=decimalLongitude",
      "Y_POSSIBLE_NAMES=decimalLatitude"
    ),
    quiet = TRUE
  )

Tanto provincias como felidae son objetos de la clase sf (y además de data.frame).

# Clase del objeto provincias
class(provincias)
[1] "sf"         "data.frame"
# Clase del objeto felidae
class(felidae)
[1] "sf"         "data.frame"

Al escribirse el nombre de un objeto sf en la consola de R, se despliega información general sobre este.

# Información general sobre el objeto provincias
provincias
Simple feature collection with 7 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 156152 ymin: 608833.8 xmax: 658879.5 ymax: 1241118
Projected CRS: CR05 / CRTM05
                 gml_id cod_catalo cod_provin  provincia
1 limiteprovincial_5k.1     160103          6 Puntarenas
2 limiteprovincial_5k.2     160103          1   San José
3 limiteprovincial_5k.3     160103          7      Limón
4 limiteprovincial_5k.4     160103          3    Cartago
5 limiteprovincial_5k.5     160103          2   Alajuela
6 limiteprovincial_5k.6     160103          5 Guanacaste
7 limiteprovincial_5k.7     160103          4    Heredia
                                                                                                                                                                                                                 ori_toponi
1 En documento de 1720, se menciona la llegada del pirata Chipperton a la zona, en el cual aparece la descripcíon referente a una embarcación pequeña en la Punta de Arena, adoptando con el tiempo el nombre de Puntarenas
2                                                                                                                                              Se remonta a la creación de la ermita dedicada al Patriarca San José en 1737
3                                                                                         El origen del nombre de la provincia se remonta a 1852, cuando por primera vez se cita en un documento oficial el puerto de Limón
4                                         Don Juan Vázques de Coronado escogió el sitio en el valle del Guarco para trasladar a la ciudad de Garcimuños, en 1563, bautizando al nuevo asentamiento con el nombre de Cartago
5                                                                                                                  Se remonta al paraje llamado La Lajuela que por primera vez se cita en los Protocolos de Cartago de 1657
6        En alegoria a un frondoso árbol de Guanacaste ubicado en la intersección de los caminos que se dirigían a Nicoya, Bagaces y Rivas, en lo que hoy día es el parque de Liberia. Esta referencia data del siglo XVIII
7                                        En correspondiencia al Presidente  de la Real Audiencia de Guatemala, Capitán General don Alonso Fernández de Heredia, de la Inmaculada Concepción de Cubujuquí a Villa de Heredia
      area  version                       geometry
1 11298.51 20201222 MULTIPOLYGON (((159917.7 60...
2  4969.73 20201222 MULTIPOLYGON (((507877.9 11...
3  9176.96 20201222 MULTIPOLYGON (((535996.1 11...
4  3093.23 20201222 MULTIPOLYGON (((520206.6 11...
5  9772.27 20201222 MULTIPOLYGON (((473981.3 11...
6 10196.32 20201222 MULTIPOLYGON (((380807.3 11...
7  2663.46 20201222 MULTIPOLYGON (((482356.4 11...

st_crs() y st_transform - manejo de sistemas de coordenadas

La función st_crs() retorna el CRS de un objeto sf.

# CRS del objeto provincias
st_crs(provincias)
Coordinate Reference System:
  User input: CR05 / CRTM05 
  wkt:
PROJCRS["CR05 / CRTM05",
    BASEGEOGCRS["CR05",
        DATUM["Costa Rica 2005",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",5365]],
    CONVERSION["Costa Rica TM 2005",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-84,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9999,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping (large and medium scale)."],
        AREA["Costa Rica - onshore and offshore east of 86°30'W."],
        BBOX[2.21,-86.5,11.77,-81.43]],
    ID["EPSG",5367]]
# CRS del objeto felidae 
st_crs(felidae)
Coordinate Reference System: NA

st_crs() también puede asignar un CRS a un objeto sf que no lo tiene.

# Asignación de un CRS al objeto felidae
st_crs(felidae) <- 4326

La función st_transform() transforma un objeto sf a un nuevo CRS.

# Transformación del CRS del objeto provincias
provincias <-
  provincias %>%
  st_transform(4326)

plot() - mapeo

La función plot() muestra objetos sf en un mapa.

# Mapeo de las geometrías del objeto provincias
plot(provincias$geometry)
# Mapeo con argumentos adicionales de plot()
plot(
  provincias$geometry,
  extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
  main = "Provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

Los argumentos reset y add de plot() permiten generar un mapa con varias capas.

# Primera capa del mapa
plot(
  provincias$geometry,
  extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
  main = "Distribución de registros de presencia de Felidae (felinos) en Costa Rica",
  axes = TRUE,
  graticule = TRUE,
  reset = FALSE
)

# Segunda capa
plot(felidae$geometry,
     add = TRUE,     
     pch = 16,
     col = "orange")

Para conocer los valores de pch, puede consultar R plot pch symbols.

st_write - escritura de datos

La función st_write() guarda en el disco un objeto sf en los diferentes formatos vectoriales de GDAL.

# Especificación del directorio de trabajo (debe utilizarse una ruta existente)
setwd("/home/mfvargas")

# Escritura del objeto provincias
provincias %>%
  st_write("provincias.shp")

# Escritura del objeto felidae
manigordos %>%
  st_write("felidae.kml")

Mapeo de objetos sf con otros paquetes

leaflet

El paquete leaflet genera mapas interactivos en HTML.

# Mapa leaflet básico con capas de provincias y registros de presencia de felinos
leaflet() %>%
  addTiles() %>%
  addPolygons(
    data = provincias,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.0,
  ) %>%
  addCircleMarkers(
    data = felidae,
    stroke = F,
    radius = 4,
    fillColor = 'orange',
    fillOpacity = 1
  )

Datos raster

El modelo raster

El modelo de datos raster usualmente consiste de un encabezado y de una matriz con celdas (también llamadas pixeles) de un mismo tamaño. El encabezado define el CRS, la extensión y el punto de origen de una capa raster. Por lo general, el origen se ubica en la esquina inferior izquierda o en la esquina superior izquierda de la matriz. La extensión se define mediante el número de filas, el número de columnas y el tamaño (resolución) de la celda.

Cada celda tiene una identificación (ID) y almacena un único valor, el cual puede ser numérico o categórico, como se muestra en la figura 2.

Figura 2. El modelo raster: (A) ID de las celdas, (B) valores de las celdas, (C) mapa raster de colores. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#raster-data)

A diferencia del modelo vectorial, el modelo raster no necesita almacenar todas las coordenadas de cada geometría (i.e. las esquinas de las celdas), debido a que la ubicación de cada celda puede calcularse a partir de la información contenida en el encabezado. Esta simplicidad, en conjunto con el álgebra de mapas, permiten que el procesamiento de datos raster sea mucho más eficiente que el procesamiento de datos vectoriales. Por otra parte, el modelo vectorial es mucho más flexible en cuanto a las posibilidades de representación de geometrías y almacenamiento de valores, por medio de múltiples elementos de datos.

Los mapas raster generalmente almacenan fenómenos continuos como elevación, precipitación, temperatura, densidad de población y datos espectrales. También es posible representar mediante raster datos discretos, tales como tipos de suelo o clases de cobertura de la tierra, como se muestra en la figura 3.

Figura 3. Ejemplos de mapas raster continuos y categóricos. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#raster-data)

El paquete terra

El paquete terra implementa un conjunto de funciones para la lectura, escritura, manipulación, análisis y modelado de datos raster.

# Lectura de una capa raster de altitud
altitud <-
  rast(
    "/vsicurl/https://raw.githubusercontent.com/tpb728O-programaciongeoespacialr/2021ii/master/datos/worldclim/altitud.tif"
  )
plot(altitud)

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY-SA 4.0. Source code is available at https://github.com/tpb728O-programaciongeoespacialr/2021ii/, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".